/*

This file makes Appendix B Figure 3

*/
global writeout 0 // change to 1 for writing out, 0 to not
global quietly quietly // comment out for noisy


$quietly { // top top quietly
  clear
  set more off
  set seed 1337
  set sortseed 1337

  capture cd "~/Desktop/RESTAT_CODE"

  global dir_root = "."
  global dir_code = "${dir_root}/code"
  global dir_data = "${dir_root}/data"
  global dir_output "${dir_root}/output"

  // From table 2 file - we used scalar variable set by that file, but here we can't assume the user will run all tables sequentially.
  scalar lambda = 16.354

  import delimited "${dir_data}\utah-history.csv", clear

  tostring date, replace
  rename date date_str
  gen date = date(date_str, "MDY")
  format date %td

  gen pos_rate = positiveincrease/totaltestresultsincrease

  rename positive ccases

  sort date, stable

  tset date
  gen pos_incr_ma7 = (positiveincrease + l1.positiveincrease + l2.positiveincrease + l3.positiveincrease + l4.positiveincrease + l5.positiveincrease + l6.positiveincrease)/7
  gen totaltest_incr_ma7 = (totaltestresultsincrease + l1.totaltestresultsincrease + l2.totaltestresultsincrease + l3.totaltestresultsincrease + l4.totaltestresultsincrease + l5.totaltestresultsincrease + l6.totaltestresultsincrease)/7
  gen pos_rate_ma7 = pos_incr_ma7 / totaltest_incr_ma7

  *ILI Data from https://www.cdc.gov/flu/weekly/weeklyarchives2020-2021/data/senAllregt11.html
  gen ili = 1
  replace ili = 1.1 if date >= mdy(09,28,2020)
  replace ili = 1.2 if date >= mdy(10,05,2020)
  replace ili = 1.1 if date >= mdy(10,12,2020)
  replace ili = 1.1 if date >= mdy(10,19,2020)
  replace ili = 1.2 if date >= mdy(11,02,2020)
  replace ili = 1.3 if date >= mdy(11,09,2020)
  replace ili = 1.4 if date >= mdy(11,16,2020)
  replace ili = 1.4 if date >= mdy(11,23,2020)
  replace ili = 1.5 if date >= mdy(11,30,2020)
  replace ili = 1.4 if date >= mdy(12,07,2020)
  replace ili = 1.4 if date >= mdy(12,14,2020)
  replace ili = 1.4 if date >= mdy(12,21,2020)
  replace ili = 1.5 if date >= mdy(12,28,2020)
  replace ili = 1.5 if date >= mdy(01,04,2021)
  replace ili = 1.4 if date >= mdy(01,11,2021)
  replace ili = 1.3 if date >= mdy(01,18,2021)
  replace ili = 1.2 if date >= mdy(01,25,2021)
  replace ili = 1.1 if date >= mdy(02,01,2021)
  replace ili = 1.1 if date >= mdy(02,07,2021)
  replace ili = 1 if date >= mdy(02,15,2021)

  *Adjust the latent prevalene estimates
  gen pI_full = (pos_rate_ma7)/(`=lambda' * (1 - pos_rate_ma7) + pos_rate_ma7) * 100000
  gen pI_full_ili = (pos_rate_ma7)/(`=lambda' / ili * (1 - pos_rate_ma7) + pos_rate_ma7) * 100000

  label variable pI_full "Latent COVID-19 Prevalence"
  label variable pI_full_ili "Latent COVID-19 Prevalence (corrected for flu season)"

  scalar sdate1 = td(1april2020)
  scalar sdate2 = td(01june2020)
  scalar sdate3 = td(01aug2020)
  scalar sdate4 = td(01oct2020)
  scalar sdate5 = td(01dec2020)
  scalar sdate6 = td(01feb2021)

  *graph new cases and our prevalence measure
    graph set window fontface "Times New Roman"

  graph twoway ///
      (line pI_full date if date > mdy(04,15,2020) & date <= mdy(02,21,2023), lwidth(med thick) lcolor(red) lpattern(dash solid)) ///
      (line pI_full_ili date if date > mdy(04,15,2020) & date <= mdy(02,21,2023), lwidth(thick) lcolor(gray) lpattern(solid)) ///
    , xlabel(`=sdate1' `=sdate2' `=sdate3' `=sdate4' `=sdate5' `=sdate6', format(%tdMon-DD-YYYY) angle(45)) ///
      xtitle("Time") ytitle("Fraction of population (per 100,000)") ///
      legend(rows(2)) ///
      ysize(4) xsize(8) scale(1.1) scheme(s2color) graphregion(color(white))  ///
      name(flue_correction, replace)

  if $writeout graph export "${dir_output}/ili_adjustment.pdf", replace name(flue_correction)

} // end top top $quietlys
